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Abstract 

In this article we show in detail the derivation of an integration scheme for the dis- 
sipative particle dynamic model (DPD) using the stochastic Trotter formula [1]. We 
explain some subtleties due to the stochastic character of the equations and exploit 
analyticity in some interesting parts of the dynamics. The DPD- Trotter integrator 
demonstrates the inexistence of spurious spatial correlations in the radial distribu- 
tion function for an ideal gas equation of state. We also compare our numerical 
integrator to other available DPD integration schemes. 
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1 Introduction 



Mesoscopic models require the use of stochastic differential equations (SDEs) 
to include the effects of thermal fluctuations so the selection of an efficient 
stochastic integration scheme is crucial to simulate correctly these systems. In 
this article we focus on dissipative particle dynamics (DPD) [2,3] as one of the 
most simple and widely used model. Although nowadays the conventional DPD 
model has a good theoretical basis, in the past few years there has been quite 
a controversy reported in literature about practical aspects of the simulations: 
the appearance of spurious effects related to time discretization, in concrete, 
the unphysical systematic drift of the temperature from the value predicted 
by the fluctuation-dissipation theorem and uncontrolled spatial correlations 
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among particles. This is the reason of the increasing interest in developing good 
integrator methods for the DPD model. Several authors [4,5,6] have considered 
improvements to the basic stochastic Euler scheme through the use of solvers 
that have been successfully employed for deterministic dynamical systems in 
molecular dynamics (MD) simulations [7] such as the velocity Verlet algorithm 
(DPD-VV [4]). 

Recently in [1] we investigated the applicability of the Trotter formula (widely 
used in molecular simulations) to a general SDEs and discuss the optimum 
way to spht the dissipative-stochastic generators. It resulted that the Trotter 
formula cannot be applied without considering the special stochastic character 
of the equations. In general, different variables depending on the same noise 
should not be split. In the DPD equations this does not happen which allowed 
us to write a integration scheme. In [1] we concluded that, considering the 
accuracy of the equilibrium temperature and the computational cost, DPD- 
Trotter is among the best integrators for the DPD equations. 

In this article we explain in detail how to apply the stochastic Trotter formula 
to the particular case of DPD. The aim is to furnish a non-trivial example 

to be used as a reference when one wishes to derive new integration schemes 
based on the stochastic Trotter formula for a general set of SDEs. We also 
test the behavior of the radial distribution function in the DPD model. For 
an ideal gas equation of state we find the DPD- Trotter scheme presents no 
spatial correlations at any scale. 



2 A Trotter integration scheme for dissipative particle dynamics 



The DPD model consists of a set of particles moving in continuous space. 
Each particle k is defined by its position r^ and its momentum p^ and mass 
m. The dynamics is specified by a set of Langevin equations very similar to 
the molecular dynamics equations, but where in addition to the conservative 
forces there are dissipative and fluctuating forces as well 



drk = ^dt, 

dPk = Hf^k^ki [{akiFcirki) - ^uj^{rki){eki ■ Pki)) dt + y/2jk^uj{rki)dWl 

(1) 

where Fc{r) is the conservative pair interaction force weighted by positive and 
symmetric parameter aki, r^i = — Vi is the distance between the particle k 
and particle /, r^i its length and 8,^./ = Vki/rki- The weight function u usually 
has a finite range Vc- A typical selection is uj{r) = 1 — r/r^ for r < Vc and 
a;(r) = for r > Tc. The conservative force is usually chosen to be of the 
form Fc{rki) — w{rki)- This system has a well defined Gibbsian equilibrium 
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state at a temperature Tq. Because the stochastic term in conventional DPD 
does not depend on the momenta, note the Ito or Stratonovich interpretation 
are exactly equivalent (additive noise) and we can apply the standard rules of 
ordinary calculus formally treating dW as "dt". 

The global state is x = (ri, ttv, Pi, Pat) and the SDEs (1) can be ex- 
pressed as (ix = C[K]dt with formal solution x(t) = Te'^*[x](xo) where T is 
the time-ordered operator (see [8]). The global time generator is divided in two 
operators generating "orthogonal dynamics" JC — jCr + ^p, where Cr — J2k 
and £p = J2k,i>k Because in the DPD model the forces between interacting 
particles k and / satisfy action-reaction (Newton's third law), the momentum 
is locally (and totally) conserved. For each partial dynamics, the generator 
can be subdivided in components = J2fj. ^'^'^ = ^^p^^ with 
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akiFc{rki) i^D{rki){^ki ■ Pki) 



S^.Jt)^f(t)^/2^oU^(r,i)e^ki(^t ' ^0' (2) 



Note that the momentum operator has two contributions, the deterministic 
and the stochastic which is an explicit function of time with f{t) = dW^i/dt. 

The formal solution of our system x(t) = Te^*[x](xo) corresponds to a con- 
tinuous time evolution. In order to devise any integrator scheme we must 
discretize the continuous time in finite steps. The continuum time propaga- 
tor can be approximated by discrete time steps of size At — t/P, recursively 

applying P times the exponential operator e^* = (e^^^^^ ~ c^"^* • • • e'^^*. 
Note that when the generator depends exphcitly on time jC,{t) = £*, the time- 
ordered exponential is relevant and the recursively nested exponentials become 
Te^** ^ g£*+ *At...g/;*+ *Atg£*+ *At Qg]), At this point we must provide 

some approximation of the discrete time propagator of a "generic" global dy- 
namics e^^*. As we have mentioned in the previous paragraph, in general the 
generator C is formed by many generators J2 J^i each of them corresponding 
to a particular dynamics i. The generalized Trotter formula (Strang [9]) gen- 
erates a straightforward approximation to the time propagator exact up to 
second order in time 

\i=M j=l J 

Because it can be performed in many possible ways, the most important prac- 
tical issue to apply formula (3) is the selection of a particular splitting of the 
global dynamics. One reasonable criteria is to keep the minimum number of 
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generators and exploit analyticity for each of them whenever possible. For the 
stochastic equations of DPD, the splitting we propose consist in 1 + 
operators: the global Cr and a Cp for each pair. 

The Baker-Campbell-Haussdorff (BCH) formula reads 

e-^e^ = e'^+^+5[-^'^]+^[-4>[AB]]+^[B,[B,^]]+- 

so for [A, <B] = we have the exact formula e'^+^ = e-^e^ = e^e"^. The DPD 
position generator Cj. — J2k u is composed by many simple individual 

generators per particle and components that satisfy [Cy.i^,Cr'^] = for all 
particles k,l {k ^ I) and components /x, u except \x v- Therefore we can use 
the exact formula 

g£.At ^ gE.^.^At ^ fj /"jj e'^^^A (5) 

fe=l \\x=\ ) 

with d the dimensionality. In MD (DPD without dissipative and random 
forces) the momentum generator Cp = J2k,fi'^p'^^ ^.Iso satisfies [>CpM ,j0.p<j^J = 
for all pairs of particles kp, Iq, with ky^l, py^q, p^q and components 
IJ>,iy, /I 7^ ly, such that the ordering of the individual-component momentum 

generators is absolutely irrelevant. On the contrary in DPD, the forces de- 
pend on the other components of the velocity of the particle and also on other 

/* Af \ ^ r.^^ /\f. 

particles velocities and the operator e = e/^k,i>k p cannot be globally 
integrated and has to be approximated in some way. This is the reason for the 
splitting it in momentum operators. Due to this splitting and formula 

(3), the DPD scheme is finally given by the following Trotter integrator 

^it + At)^( n e^'^''^](f[e^rAt\( e^^'^)x(i). (6) 

\g=l,r>l / \i=l / \k=N,l<N J 

The propagator that corresponds to the generator produces the position 
update which is analytically given by 

e< [x]: r,-(t + At) = + (7) 

because the momentum is a constant in this step of the scheme. The next 
step is to solve the propagator of the momenta of the interaction pair /c, / 
(corresponding to the generator C^^) independently of the positions. We have 
mentioned before that DPD forces satisfy action-reaction, so for a particular 
interacting pair k.l we propose to make a change of variables from Pa;,P; 
to Pit + Pi,P/fc/ = Pfc — Pz- The new system to solve is d{'p^. + pz) = and 
d^ki = 2(ipA;. Because the positions of the particles are "frozen" at this step 
of the Trotter scheme, the equation for dpki can be solved more easily for the 
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projection on the radial direction p^i — pki ■ eu 

dpli = Adt-Bplidt + CdWli, (8) 

where A = 2akiFc{rki), B = 2'j/muj^ and C = 2y^2rfkBToUJ. This equation is 
an Ornstein-Uhlenbeck process with analytical solution [10] 

Pliit) = e-^^YM + Af e^(-*)ds + C /* e^(-*)diy„ (9) 

Jto Jto 

where At — t — to, to being the initial time. The solution of (9) requires the 
generation of colored noise based on a numerical scheme itself. A version of 
the method to generate coloured noise [1,11] adapted to Eq.(9) results 

ApL = (p.. • e.. - ^) (e-^ - l) + ^2k,Tom (l - 6"^)^^ (10) 

where r = ^jmiJ^, = are normal distributed with zero mean and 
variance one (A^(0, 1)) and Ap|; = ~ Pw(^o)- The propagator e'-'p ^* for 

Pfe and Pi gives 

e^^'^lx] : (p^^*, pr^O = (p. W + ^e.Kt), pK^) " ^e.^t)) • (H) 

So in DPD we can solve the dynamics corresponding to the generator L^^ 
(globally for all components at the same time) without the need to go to the 
scalar operator corresponding to the coordinate \i. In practice the Trotter 
integration algorithm (6) consists of the following steps: for the interaction 
pairs I update the momentum half timestep according to the propagator 
(11) with a noise iterate over particles k updating the position according 
to (7); finally, update pairs k,l in reverse order again using the propagator 
(11) but with new noises This algorithm requires the calculation of the 
pair-list only once per iteration and has the same complexity as a simple DPD 
velocity- Verlet scheme (DPD-VV [4]). 

We tested in [1] this integration scheme using the open-source code mydpd 
[12] for the equilibrium temperature with = 4000 particles, 7 = 4.5, /c^To = 
1, m = 1, Tc = 1 in a three dimensional periodic box {L,L,L) with L = 10 
with periodic boundary conditions. These settings give a particle density p — 
4. Here, we show in left Fig.l the radial distribution function for an = 
(corresponding to an ideal gas equation of state) and a time step At = 0.05. 
We compare the results for three methods: the velocity Verlet (DPD-VV) [4], 
the Shardlow scheme [13] and DPD- Trotter [1] . We find good agreement with 
the theoretical value 1 for Shardlow and DPD- Trotter integrators but DPD- 
VV is notably wrong displaying spurious spatial correlations at distances less 
than the finite range Vf.. In right Fig. 1 we show the radial distribution function 
for a simulation with a^i = 25 and a time step At = 0.01. As we see the three 
methods perform very similarly. 
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Fig. 1. Radial distribution function for three integrator methods. Velocity Verlet 
with (A) symbols, Shadlow scheme with (□) symbols and Trotter DPD with (•) 
symbols. Left figure corresponds to an ideal gas simulation aki = 0. Right figure 
corresponds to simulations including conservative forces with Uki = 25. 

3 Conclusions 



The stochastic Trotter formula can be successfully applied to the DPD model 
and the procedure to tailor the integrator scheme has been explained in detail. 
In the scheme we have also exploited the exact integration of important parts 
of the dynamics like the conservation of total momentum of an interacting 
pair of particles. The DPD- Trotter integrator displays correctly the radial 
distribution functions for an ideal gas (no conservative forces among particles) 
and also for a non ideal gas. Following this important example and [1] it should 
be strathforward to apply the stochastic Trotter formula to new mesoscopic 
models and more general SDEs. 
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